
#define kmax 19
#define ksqmax 250
#define nmax 2
#define nsqmax 16
#define maxk 10000

class ewald
{
	/******************************************************************
    ** REAL-SPACE AND RECIPROCAL-SPACE PARTS OF EWALD SUM FOR IONS.  **
    **                                                               **
    ** REFERENCES:                                                   **
    **                                                               **
    ** WOODCOCK AND SINGER, TRANS. FARADAY SOC. 67, 12, 1971.        **
    ** DE LEEUW ET AL., PROC. ROY. SOC. A 373, 27, 1980.             **
    ** HEYES, J. CHEM. PHYS. 74, 1924, 1981.                         **
    ** SEE ALSO FINCHAM, MDIONS, CCP5 PROGRAM LIBRARY.               **
    **                                                               **
    ** ROUTINES SUPPLIED:                                            **
    **                                                               **
    ** setup()														 **
    **    SETS UP THE WAVEVECTORS FOR USE IN THE EWALD SUM           **
    ** r_ewald ()					                                 **
    **    CALCULATES THE R-SPACE PART OF THE SUM                     **
    ** k_ewald ()					                                 **
    **    CALCULATES THE K-SPACE PART OF THE SUM                     **
    ** erfc(x)				                                         **
    **    RETURNS THE COMPLEMENTARY ERROR FUNCTION                   **
    **                                                               **
    ** PRINCIPAL VARIABLES:                                          **
    **                                                               **
    ** int  totk	         THE TOTAL NUMBER OF K-VECTORS STORED    **
    ** int  maxk	         MAXIMUM POSSIBLE NUMBER OF K-VECTORS    **
    ** int  kmax	         MAX INTEGER COMPONENT OF THE K-VECTOR   **
    ** int  ksqmax	         MAX SQUARE MOD OF THE K-VECTOR REQUIRED **
    ** double vr             ENERGY FROM R-SPACE SUM                 **
    ** double vk	         ENERGY FROM K-SPACE SUM                 **
    ** double kvec[maxk]	 ARRAY USED TO STORE K-VECTORS           **
    ** double   kappa        WIDTH OF CANCELLING DISTRIBUTION        **
    **                                                               **
    ** USAGE:                                                        **
    **                                                               **
    ** setup IS CALLED ONCE AT THE BEGINNING OF THE SIMULATION       **
    ** TO CALCULATE ALL THE K-VECTORS REQUIRED IN THE EWALD SUM.     **
    ** THESE VECTORS ARE USED THROUGHOUT THE SIMULATION IN THE       **
    ** SUBROUTIN k_ewaldTO CALCULATE THE K-SPACE CONTRIBUTION TO THE **
    ** POTENTIAL ENERGY AT EACH CONFIGURATION. THE SELF TERM IS      **
    ** SUBTRACTED FROM THE K-SPACE CONTRIBUTION IN k_ewald.          **
    ** THE SURFACE TERM FOR SIMULATIONS IN VACUUM IS NOT INCLUDED.   **
    ** ROUTINE r_ewald RETURNS THE R-SPACE CONTRIBUTION TO THE EWALD **
    ** SUM AND IS CALLED FOR EACH CONFIGURATION IN THE SIMULATION.   **
    ** A CUBIC BOX AND UNIT BOX LENGTH ARE ASSUMED THROUGHOUT.       **
    *******************************************************************/
private:
	double twopi;
	double twopioverbox;
public:

	double kappa;
	double *kvec;
	double B;
	int totk;
	//int npart;
	double **pos;
	double vr,vk;
	double phi;
	double lbox;
	double rsqpi;
	double volume;
	double selfk;

	ewald(double kap,double lb,  double **posiciones)
	{
		pos=posiciones;
		//npart=np;
		kappa=kap/lb;
		twopi=2.0*pi;
		rsqpi = 0.5641896;//1./sqrt(pi)
		lbox=lb;
		volume=lb*lb*lb;
		twopioverbox=twopi/lbox;
		B=1.0 / 4.0 / kappa / kappa;

		kvec = new double [maxk];
		selfk=2.*pi/3.*pow(kappa/pi,3./2.)*npart;
	}

	~ewald()
	{
		delete []kvec;
		//for (int i=0;i<npart;i++)
		//{
		//	delete pos[i];
		//}
		//delete []pos;
	}

	void setup(void);
	void setupdip(void);
	void r_ewald(void);
	void r_ewaldvar(void);
	void k_ewald(void);
	void r_ewaldDip(void);
	void k_ewaldDip(void);
	double distancia2(int i,int j);
	double distancia(int i,int j);
	double distancia2(int i,int j, double *dz);
	double distancia(int i,int j,double *dz);
	double distancia2Comp(int ind1,int ind2, double *dz,int nx,int ny, int nz);
	double distanciaComp(int ind1,int ind2,double *dz, int nx,int ny, int nz);
	double erfc(double x);
	double edirecta(void);
	double edirectaDipolar(void);
	void BCFuntions(double r, double *B, double *C);
	double edirectaDipolarReplicas();
};

void ewald::setup(void)
{
   /* *******************************************************************
    ** ROUTINE TO SET UP THE WAVE-VECTORS FOR THE EWALD SUM.         **
    **                                                               **
    ** THE WAVEVECTORS MUST FIT INTO A BOX OF UNIT LENGTH.           **
    ** IN THIS EXAMPLE WE ALLOW A MAXIMUM OF 1000 WAVEVECTORS.       **
    ********************************************************************/

	double rkx, rky, rkz, ksq, rksq;
	totk=0;
	//** LOOP OVER K-VECTORS. NOTE KX IS NON-NEGATIVE **

	for(int kx=0;kx<=kmax;kx++)
	{
		rkx=twopioverbox*double(kx);
		for(int ky=-kmax;ky<=kmax;ky++)
		{
			rky=twopioverbox*double(ky);

			for(int kz=-kmax;kz<=kmax;kz++)
			{
				rkz=twopioverbox*double(kz);
				ksq=kx*kx+ky*ky+kz*kz;

				if(ksq<ksqmax && ksq!=0)
				{

					if ( totk>maxk)
					{
						cout<<"kvec es demasiado pequeño"<<endl;
						exit(1);
					}
					rksq=rkx*rkx+rky*rky+rkz*rkz;
					kvec[totk]=twopi*exp(-B*rksq)/rksq/volume;
					//cout<<totk<<"  "<<ksq<<"  "<<kx<<"  "<<ky<<"  "<<kz<<"  "<<kvec[totk]<<endl;
					totk++;
				}

			}
		}
	}
    //cout<<"Suma de Ewald completa"<<endl;
	//cout<<"El numero de vectores de onda considerado es de "<<totk<<endl;
}

void ewald::setupdip(void)
{
   /* *******************************************************************
    ** ROUTINE TO SET UP THE WAVE-VECTORS FOR THE EWALD SUM.         **
    **                                                               **
    ** THE WAVEVECTORS MUST FIT INTO A BOX OF UNIT LENGTH.           **
    ** IN THIS EXAMPLE WE ALLOW A MAXIMUM OF 1000 WAVEVECTORS.       **
    ********************************************************************/

	double rkx, rky, rkz, ksq, rksq;
	totk=0;
	//** LOOP OVER K-VECTORS. NOTE KX IS NON-NEGATIVE **

	for(int kx=0;kx<kmax;kx++)
	{
		rkx=twopioverbox*double(kx);
		for(int ky=-kmax;ky<kmax;ky++)
		{
			rky=twopioverbox*double(ky);

			for(int kz=-kmax;kz<kmax;kz++)
			{
				rkz=twopioverbox*double(kz);
				ksq=kx*kx+ky*ky+kz*kz;

				if(ksq<ksqmax && ksq!=0)
				{
					//totk++;
					if ( totk>maxk)
					{
						cout<<"kvec es demasiado pequeño"<<endl;
						exit(1);
					}
					rksq=rkx*rkx+rky*rky+rkz*rkz;
					//kvec[totk]=twopi*exp(-B*rksq)/ksq/volume;
					kvec[totk]=twopi*exp(-B*rksq)/rksq/volume ;
					totk++;
				}

			}
		}
	}
    //cout<<"Suma de Ewald completa"<<endl;
	//cout<<"El numero de vectores de onda considerado es de "<<totk<<endl;
}

double ewald::distancia2(int ind1,int ind2)
{
	double dr, dis2, disaux;
	double an;
	dis2=0.;
	disaux=0.;
	for(int k=0;k<3;k++)
	{
		dr=pos[ind2][k]-pos[ind1][k];
		disaux+=dr*dr;
		dr-=lbox*anint(dr/(lbox));
		dis2+=dr*dr;
	}
        return dis2;
}

double ewald::distancia(int ind1,int ind2, double *dz)
{
	double xdz;
	double d2=distancia2(ind1,ind2,&xdz);
	*dz=xdz;
	return sqrt(d2);
}

double ewald::distancia2(int ind1,int ind2, double *dz)
{
	double dr, dis2;
	dis2=0.;
	for(int k=0;k<3;k++)
	{
		dr=pos[ind1][k]-pos[ind2][k];
		dr-=lbox*anint(dr/(lbox));
		dis2+=dr*dr;
	}
	*dz=dr;
    return dis2;
}

double ewald::distancia2Comp(int ind1,int ind2,double *dz, int nx,int ny, int nz)
{
	double dis2;
	double dr[3];
	int dn[3]={nx,ny,nz};
	dis2=0.;
	for(int k=0;k<3;k++)
	{
		dr[k]=pos[ind2][k]-pos[ind1][k]+dn[k]*lbox;
		//dr-=lbox*anint(dr/(lbox));
		dis2+=dr[k]*dr[k];
	}
	*dz=dr[2];
    return dis2;
}

double ewald::distancia(int ind1,int ind2)
{
	double d2=distancia2(ind1,ind2);
	return sqrt(d2);
}

double ewald::distanciaComp(int ind1,int ind2, double *dz,int nx,int ny, int nz)
{
	double d2=distancia2Comp(ind1,ind2,dz,nx,ny,nz);
	return sqrt(d2);
}

void ewald::r_ewald(void)
{
    /******************************************************************
    ** CALCULATES R-SPACE PART OF POTENTIAL ENERGY BY EWALD METHOD.  **
    *******************************************************************/
	double dr[3], rij;
	double z1=1.0;
	double z2=1.0;
	double krij,vij;
	vr=0.0;

	for(int i=0;i<npart-1;i++)
	{
		for(int j=i+1;j<npart;j++)
		{
			rij=distancia(i,j);
			krij=kappa*rij;
			vij = erfc ( krij ) / rij;
			vr+=vij;
		}
	}
}

void ewald::r_ewaldDip(void)
{
    /******************************************************************
    ** CALCULATES R-SPACE PART OF POTENTIAL ENERGY BY EWALD METHOD.  **
    *******************************************************************/
	double dr[3], rij;
	double krij,vij, dz;
	double B,C;
	vr=0.0;

	for(int i=0;i<npart-1;i++)
	{
		for(int j=i+1;j<npart;j++)
		{
//			rij=distancia(i,j,&dz);
//			BCFuntions(rij, &B, &C);
//			vij=B-dz*dz*C;
//			vr+=vij;

			rij=distanciaComp(i,j,&dz,0,0,0);
			BCFuntions(rij, &B, &C);
			vij=B-dz*dz*C;
			vr+=vij;

			rij=distanciaComp(i,j,&dz,1,0,0);
			BCFuntions(rij, &B, &C);
			vij=B-dz*dz*C;
			vr+=vij;

			rij=distanciaComp(i,j,&dz,-1,0,0);
			BCFuntions(rij, &B, &C);
			vij=B-dz*dz*C;
			vr+=vij;

			rij=distanciaComp(i,j,&dz,0,1,0);
			BCFuntions(rij, &B, &C);
			vij=B-dz*dz*C;
			vr+=vij;

			rij=distanciaComp(i,j,&dz,0,-1,0);
			BCFuntions(rij, &B, &C);
			vij=B-dz*dz*C;
			vr+=vij;

			rij=distanciaComp(i,j,&dz,0,0,1);
			BCFuntions(rij, &B, &C);
			vij=B-dz*dz*C;
			vr+=vij;

			rij=distanciaComp(i,j,&dz,0,0,-1);
			BCFuntions(rij, &B, &C);
			vij=B-dz*dz*C;
			vr+=vij;
		}
	}
}

void ewald::BCFuntions(double r, double *B, double *C)
{
	double ar = kappa*r;
	double ar2= ar*ar;
	double arpi=2.0*ar/sqrt(pi);
	double r3=r*r*r;
	double r5=r3*r*r;
	double auxB, auxC;

	double erc=erfc(ar);
	auxB=(erc+arpi*exp(-ar2))/r3;
	auxC=(3.0*erc+arpi*(3.+2.*ar2)*exp(-ar2))/r5;
	*B=auxB;
	*C=auxC;
}

void ewald::k_ewald(void)
{
    /******************************************************************
    ** CALCULATES K-SPACE PART OF POTENTIAL ENERGY BY EWALD METHOD.  **
    **                                                               **
    ** THE SELF TERM IS SUBTRACTED.                                  **
    ** IN ONE COORDINATE DIRECTION (X), SYMMETRY IS USED TO REDUCE   **
    ** THE SUM TO INCLUDE ONLY POSITIVE K-VECTORS.                   **
    ** THE NEGATIVE VECTORS IN THIS DIRECTION ARE INCLUDED BY USE    **
    ** OF THE MULTIPLICATIVE VARIABLE 'FACTOR'.                      **
    **                                                               **
    ** PRINCIPAL VARIABLES:                                          **
    **                                                               **
    ** INTEGER N                   NUMBER OF IONS                    **
    ** REAL    RX(N),RY(N),RZ(N)   POSITIONS OF IONS                 **
    ** REAL    Z(N)                IONIC CHARGES                     **
    ** REAL    VK                  K-SPACE POTENTIAL ENERGY          **
    ** REAL    VKS                 SELF PART OF K-SPACE SUM          **
    ******************************************************************/

	complex <double>sum;
	double factor, ksq;
	double vd,vs;
	complex<double> eikx[npart][kmax+1];
	complex<double> eiky[npart][2*kmax+1];
	complex<double> eikz[npart][2*kmax+1];
	complex<double> eikr[npart];

	//CONSTRUCT EXP(IK.R) FOR ALL IONS AND K-VECTORS **
	//CALCULATE KX, KY, KZ = 0 , -1 AND 1 EXPLICITLY **

	for(int i=0;i<npart;i++)
	{
		eikx[i][0]=complex<double>(1.0,0.0);
		eiky[i][kmax+0] = complex<double>(1.0, 0.0);
        eikz[i][kmax+0] = complex<double>(1.0, 0.0);
		eikx[i][1]      = complex<double>(cos(twopioverbox*pos[i][0]) ,sin(twopioverbox*pos[i][0]));
		eiky[i][kmax+1] = complex<double>(cos(twopioverbox*pos[i][1]) ,sin(twopioverbox*pos[i][1]));
        eikz[i][kmax+1] = complex<double>(cos(twopioverbox*pos[i][2]) ,sin(twopioverbox*pos[i][2]));
		eiky[i][kmax-1]= conj(eiky[i][kmax+1]);
        eikz[i][kmax-1]= conj(eikz[i][kmax+1]);
	}

	//** CALCULATE REMAINING KX, KY AND KZ BY RECURRENCE **

	for (int kx=2;kx<=kmax;kx++)
	{
		for(int i=0;i<npart;i++)
		{
			eikx[i][kx]=eikx[i][kx-1]*eikx[i][1];
		}
	}

	for(int ky=2;ky<=kmax;ky++)
	{
		for(int i=0;i<npart;i++)
		{
			eiky[i][kmax+ky]=eiky[i][kmax+ky-1]*eiky[i][kmax+1];
			eiky[i][kmax-ky]=conj(eiky[i][kmax+ky]);

			eikz[i][kmax+ky]=eikz[i][kmax+ky-1]*eikz[i][kmax+1];
			eikz[i][kmax-ky]=conj(eikz[i][kmax+ky]);
		}
	}

	ofstream MiArchivo ("k:\\reciproco.dat",ios::out);
	MiArchivo.is_open();
	//** SUM OVER ALL VECTORS **
	vd=0.0;
	totk=0;

	for(int kx=0;kx<=kmax;kx++)
	{
		if(kx==0)
		{
			factor=1.0;
		}
		else
		{
			factor=2.0;
		}
		for(int ky=-kmax;ky<=kmax;ky++)
		{
			for(int kz=-kmax;kz<=kmax;kz++)
			{
				ksq =kx*kx+ky*ky+kz*kz;
				if((ksq<ksqmax)&&(ksq!=0))
				{
					sum= complex<double>(0.0, 0.0);
					for(int i=0;i<npart;i++)
					{
						 eikr[i]=eikx[i][kx]*eiky[i][ky+kmax]*eikz[i][kz+kmax];
						 sum+=eikr[i];
					}
					vd +=factor*real(kvec[totk]*conj(sum)*sum);
					MiArchivo<<kx<<"  "<<ky<<"  "<<kz<<"  "<<factor*real(kvec[totk]*conj(sum)*sum)<<"  "<<vd<<endl;

					//cout<<totk<<"  "<<vd<<endl;
					totk++;
				}
			}
		}
	}
	MiArchivo.close();
	//** CALCULATES SELF PART OF K-SPACE SUM **
	vs=rsqpi*kappa*npart;
	//vs=0.;
	//** CALCULATE THE TOTAL K-SPACE POTENTIAL **
	vk=vd-vs;
}

void ewald::k_ewaldDip(void)
{
    /******************************************************************
    ** CALCULATES K-SPACE PART OF POTENTIAL ENERGY BY EWALD METHOD.  **
    **                                                               **
    ** THE SELF TERM IS SUBTRACTED.                                  **
    ** IN ONE COORDINATE DIRECTION (X), SYMMETRY IS USED TO REDUCE   **
    ** THE SUM TO INCLUDE ONLY POSITIVE K-VECTORS.                   **
    ** THE NEGATIVE VECTORS IN THIS DIRECTION ARE INCLUDED BY USE    **
    ** OF THE MULTIPLICATIVE VARIABLE 'FACTOR'.                      **
    **                                                               **
    ** PRINCIPAL VARIABLES:                                          **
    **                                                               **
    ** INTEGER N                   NUMBER OF IONS                    **
    ** REAL    RX(N),RY(N),RZ(N)   POSITIONS OF IONS                 **
    ** REAL    Z(N)                IONIC CHARGES                     **
    ** REAL    VK                  K-SPACE POTENTIAL ENERGY          **
    ** REAL    VKS                 SELF PART OF K-SPACE SUM          **
    ******************************************************************/

	complex <double>sum;
	double factor, ksq;
	double vd,vs;
	complex<double> eikx[npart][kmax+1];
	complex<double> eiky[npart][2*kmax+1];
	complex<double> eikz[npart][2*kmax+1];
	complex<double> eikr[npart];

	//CONSTRUCT EXP(IK.R) FOR ALL IONS AND K-VECTORS **
	//CALCULATE KX, KY, KZ = 0 , -1 AND 1 EXPLICITLY **

	for(int i=0;i<npart;i++)
	{
		eikx[i][0]=complex<double>(1.0,0.0);
		eiky[i][kmax+0] = complex<double>(1.0, 0.0);
        eikz[i][kmax+0] = complex<double>(1.0, 0.0);
		eikx[i][1]      = complex<double>(cos(twopioverbox*pos[i][0]) ,sin(twopioverbox*pos[i][0]));
		eiky[i][kmax+1] = complex<double>(cos(twopioverbox*pos[i][1]) ,sin(twopioverbox*pos[i][1]));
        eikz[i][kmax+1] = complex<double>(cos(twopioverbox*pos[i][2]) ,sin(twopioverbox*pos[i][2]));
		eiky[i][kmax-1]= conj(eiky[i][kmax+1]);
        eikz[i][kmax-1]= conj(eikz[i][kmax+1]);
	}

	//** CALCULATE REMAINING KX, KY AND KZ BY RECURRENCE **

	for (int kx=2;kx<=kmax;kx++)
	{
		for(int i=0;i<npart;i++)
		{
			eikx[i][kx]=eikx[i][kx-1]*eikx[i][1];
		}
	}

	for(int ky=2;ky<=kmax;ky++)
	{
		for(int i=0;i<npart;i++)
		{
			eiky[i][kmax+ky]=eiky[i][kmax+ky-1]*eiky[i][kmax+1];
			eiky[i][kmax-ky]=conj(eiky[i][kmax+ky]);

			eikz[i][kmax+ky]=eikz[i][kmax+ky-1]*eikz[i][kmax+1];
			eikz[i][kmax-ky]=conj(eikz[i][kmax+ky]);
		}
	}

	//** SUM OVER ALL VECTORS **
	vd=0.0;
	totk=0;

	for(int kx=0;kx<=kmax;kx++)
	{
		if(kx==0)
		{
			factor=1.0;
		}
		else
		{
			factor=2.0;
		}
		for(int ky=-kmax;ky<=kmax;ky++)
		{
			for(int kz=-kmax;kz<=kmax;kz++)
			{
				ksq =kx*kx+ky*ky+kz*kz;
				if((ksq<ksqmax)&&(ksq!=0))
				{
					sum= complex<double>(0.0, 0.0);
					for(int i=0;i<npart;i++)
					{
						 eikr[i]=eikx[i][kx]*eiky[i][ky+kmax]*eikz[i][kz+kmax];
						 sum+=eikr[i];
					}
					//cout<<ksq<<"  "<<sum<<endl;
					//double val = factor*kz*kz*real(kvec[totk]*conj(sum)*sum);
					double val = factor*real(kvec[totk]*conj(sum)*sum);
                    vd +=val;
					totk++;
				}
			}
		}
	}

	//** CALCULATES SELF PART OF K-SPACE SUM **
	vs=-selfk;//*npart;
	//vs=0.;
	//** CALCULATE THE TOTAL K-SPACE POTENTIAL **
	vk=vd-vs;
}

double ewald::erfc(double x)
{
    /******************************************************************
    ** APPROXIMATION TO THE COMPLEMENTARY ERROR FUNCTION             **
    **                                                               **
    ** REFERENCE:                                                    **
    **                                                               **
    ** ABRAMOWITZ AND STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS,    **
    **    NATIONAL BUREAU OF STANDARDS, FORMULA 7.1.26               **
    ******************************************************************/

	double A1 = 0.254829592;
	double A2 = -0.284496736;
	double A3 = 1.421413741;
	double A4 = -1.453152027;
	double A5 = 1.061405429;
	double P  =  0.3275911;
	double t, xsq, tp;

	t= 1.0/(1.0+P*x);
	xsq=x*x;
	tp= t*(A1+t*(A2+t*(A3+t*(A4+t*A5))));
	return tp*exp(-xsq);
}

double ewald::edirecta(void)
{
	double dist, en;
	en=0.;
	for(int i=0;i<npart-1;i++)
	{
		for(int j=i+1;j<npart;j++)
		{
			dist = distancia(i,j);
			en+=1.0/dist;
		}
	}
	return en;
}

double ewald::edirectaDipolar(void)
{
	double dist, en, dz;
	en=0.;
	for(int i=0;i<npart-1;i++)
	{
		for(int j=i+1;j<npart;j++)
		{
			dist = distancia(i,j,&dz);
			en+=1.0/(dist*dist*dist)*(1.-3.*(dz/dist)*(dz/dist));
		}
	}
	return en;
}


double ewald::edirectaDipolarReplicas(void)
{
    /**********************
	** compute sum
	*********************/
	double E;
	double qij;
	double nq;
	int nx,ny,nz;
	double r,r2,rx,ry,rz,rnx,rny,rnz;
	double Co;
	int NCUT=30;
	double NCUT2=NCUT*NCUT;
	double mad=0.;
	double coef;
	E=0.0;


	for(int i=0; i<npart; i++){
		for(int j=0; j<npart; j++){
			rx=pos[i][0]-pos[j][0];
			ry=pos[i][1]-pos[j][1];
			rz=pos[i][2]-pos[j][2];

			for(nx=-NCUT;nx<=NCUT;nx++){
				for(ny=-NCUT;ny<=NCUT;ny++){
					for(nz=-NCUT;nz<=NCUT;nz++){
						nq=nx*nx+ny*ny+nz*nz;
						if(nq < NCUT2){
							if((nq!=0) || (nq==0 &&i!=j))
							{
                                rnx=rx+nx*lbox;
                                rny=ry+ny*lbox;
                                rnz=rz+nz*lbox;
                                r2=rnx*rnx+rny*rny+rnz*rnz;
                                r=sqrt(r2);
                                Co=rnz/r;
                                E += (1.-3.0*Co*Co)/(r2*r);
    //							force[i][0]+=coef*rnx;
    //							force[i][1]+=coef*rny;
    //							force[i][2]+=coef*rnz;
							}
						}
					}
				}
			}
		}
	}

	E *= 0.5;

	return E;
}
